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High-Order  CESE  Methods  for  the  Euler  Equations 


David  L.  Bilyeu,*  Yung- Yu  Chen,  t  and  S.-T.  John  Yu  * 

The  Ohio  State  University,  Columbus,  OH  43210,  USA 


Recently,  Chang  reported  a  new  class  of  high-order  CESE  methods  for  solving  nonlin¬ 
ear  hyperbolic  partial  differential  equations.  A  series  of  high-order  algorithms  have  been 
developed  based  on  a  systematic,  recursive  formulation  that  achieves  fourth-,  sixth-,  and 
eighth-order  accuracy.  The  new  high-order  CESE  method  shares  many  favorable  attributes 
of  the  original  second-  order  CESE  method,  including:  (i)  compact  mesh  stencil  involving 
only  the  immediate  mesh  nodes  surrounding  the  node  where  the  solution  is  sought,  (ii)  the 
CFL  stability  constraint  remains  to  be  the  same,  i.e.,  <  1,  as  compared  to  the  original 
second-order  method,  and  (iii)  superb  shock  capturing  capability  without  using  an  approx¬ 
imate  Riemann  solver.  The  new  algorithm  has  been  demonstrated  by  solving  Burger’s 
equation. 

In  the  present  paper,  we  extend  Chang’s  high-order  method  for  system  of  linear  and 
nonlinear  hyperbolic  partial  differential  equations.  A  general  formulation  is  presented  for 
solving  the  coupled  equations  with  arbitrarily  high-order  accuracy.  To  demonstrate  the 
formulation,  several  linear  and  nonlinear  cases  are  reported.  First,  we  solve  a  convection 
equation  with  source  term  and  the  linear  acoustics  equations.  We  then  solve  the  Euler 
equations  for  acoustic  waves,  a  blast  wave,  and  Shu  and  Osher’s  test  case  for  acoustic  waves 
interacting  with  a  shock.  Numerical  results  show  higher-order  convergence  by  continuous 
mesh  refinement. 


I.  Introduction 

In  this  work,  we  extend  Chang’s  fourth-order  CESE  method  for  one  nonlinear  hyperbolic  equation  to  a 
system  of  coupled  hyperbolic  partial  differential  equations  (PDEs).  The  new  formulation  is  general  and  can 
be  used  to  achieve  arbitrarily  order  of  convergence.  To  demonstrate  the  capabilities  of  the  new  scheme,  we 
apply  the  method  to  solve  three  sets  of  equations:  (i)  the  one-dimensional  Euler  equations,  (ii)  the  linearized 
acoustic  equations,  and  (iii)  a  convection  equation  with  a  source  term. 

The  original  second-order  CESE  method  solves  the  hyperbolic  PDEs  by  discretizing  the  space-time 
domain  by  using  the  conservation  elements  (CEs)  and  solution  elements  (SEs).  The  profiles  of  unknowns  are 
prescribed  by  assumed  discretization  inside  SEs.  Aided  by  the  approximation  for  the  unknowns  in  the  SEs, 
space-time  flux  conservation  can  be  enforced  over  each  CE.  The  calculation  of  space-time  flux  conservation 
results  in  the  formulation  for  updating  the  unknowns  in  the  time  marching  process.  The  special  features 
of  the  CESE  method  include:  (i)  The  a  scheme,  the  core  scheme  of  the  CESE  method,  is  non-dissipative. 
(ii)  The  CESE  method  has  the  most  compact  mesh  stencil  and  it  involves  only  the  immediate  neighboring 
mesh  points  that  surround  the  mesh  node  where  the  solution  is  sought,  (iii)  The  method  uses  explicit 
integration  in  time  marching.  The  stability  criterion  is  CFL  <  1.  (iv)  No  approximate  Riemann  solver  is 
used  and  the  scheme  is  simple  and  efficient. 

This  paper  is  organized  as  the  following.  Section  II  reports  the  fourth-order  CESE  method  for  the 
coupled  equations  formulated  in  a  vector  form.  Section  III  shows  the  application  of  the  general  formulation 
to  the  one-dimensional  Euler,  linear  acoustic,  and  convection  equation.  Section  IV  provides  the  results  and 
discussions.  In  Section  V,  we  draw  our  conclusions. 
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+ Associate  Professor,  Dept,  of  Mechanical  Engineering,  yu.274@osu.edu  and  AIAA  Member. 


1  of  14 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


American  Institute  of  Aeronautics  and  Astronautics 


II.  Arbitrary-Order,  One-Dimensional  CESE  Method 


Consider  a  system  of  coupled  convection  equations: 


<9U  dF  _ 
dt  dx 


(1) 


where 


U  =f  (ui,u2,u3,  ■  ■  ■  ,um)T , 

F  =  (/i,/2,/3,  -  -  -  ,/m)T, 

S  =f  (si,  s2,  s3,  •  •  •  ,sm)T, 

There  are  M  equations  to  be  solved  in  the  system  Eq.  (1). 

The  space-time  stencil  used  in  this  derivation  is  the  same  as  that  reported  by  Chang  and  is  repeated 
here  for  completeness.  In  Fig.  (1)  the  solid  dots  A,  C ,  and  E  are  the  solution  points.  A  is  at  (xj,tn)  and  C 
and  E  are  at  (xj_1/2,  tn~1^2)  and  (xj+1/2,  tn_1^2),  respectively.  P+  is  between  M+  and  F.  P~  is  between 
M~  and  B.  The  distance  between  P±  and  M±  is  determined  by  a  parameter  r.  The  rectangles  ABCD  and 
ADEF  are  basic  CEs  (BCEs),  while  the  rectangle  BCEF  is  the  compound  CE  (CCE)  associated  with  the 
solution  point  A. 


B  P~  M~  A  M+  P+  F 


To  facilitate  the  discussion,  we  let  SE(j,  n)  denotes  the  SE  located  at  x  =  Xj  and  t  =  tn.  To  denote 
high-order  derivatives,  we  use  the  following  notations: 

da+bUm 

Umxatb  dxnb 

In  SE(j,  n),  the  unknown  variables  um,  m  =  1, . . . ,  M,  are  approximated  by  a  Taylor  series: 


Nm  Nm—cl  (  \n 

def  \  \  >  \V"mxatb)j 


a= 0 


E 

fe=0 


alb! 


(x-xj)a(t-tny 


(2) 


where  Nm  is  the  desired  order  subtracted  by  1,  e.g.,  for  the  fourth-order  scheme,  Nm  =  3.  The  superscript 
*  represents  the  numerical  approximation  of  the  variable.  Inside  of  a  SE  umxat t  are  constant.  The  flux 
functions  /m,  m  =  1, . . . ,  M,  can  also  be  represented  with  the  Taylor  expansion  as: 


Nm  NM—a 

a— 0  6—0 


(frr 


c“tb) 


alb ! 


(x-Xjr  c t-nb 


(3) 
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Inside  a  SE,  fmxa.tb  are  constant.  An  advantage  of  a  Taylor  series  is  that  its  derivatives  can  also  be  expressed 
as  a  Taylor  series 


A  B-a 


—  CL  /  \  fi 

■  (*» j. *0  =  EE  uT+t  J  ( x - xo)a (t - *n)fa ;  BA=NM-z-i 


a— 0  6=0 


a!6! 


and 


A  B-o 


(  fmxa+ztb+i ) 7 


=  E  E  ^  ~ Xj^a ^ ~ r')b ;  sIjvm" z~i- 

a=0  6=0 


(4) 


(5) 


Equations  2  and  3  contain  2  n  unknowns.  In  the  following  derivation,  it  will  be  shown,  for 

a  given  SE (j,n),  the  only  independent  variables  are  the  spatial  derivatives  (zim)",  ( umx )",  ...,  ( umX“ )", 
to  =  1, ... ,  M,  a  =  IVm-  Since  it  is  assumed  that  the  fluxes  are  known  functions  of  the  flow  variables.  The 
flux  terms  in  Eq.  (3)  can  be  determined  from  the  chain  rule.  To  proceed,  we  define: 


fmi  — 


def  dfn 


dui  duiduk 

where  to,  l,  k,p  =  1, . . . ,  M.  For  SE(j,  n),  we  obtain: 

Neq 


f  def  d2fm  f  def  <93  /„ 

Jmi  k  o  o  j  Jn 


1  J  7Yll  ,k  ,p 


duidukdup 


(6) 


dy  i 

d2fm 

dyidy2 

d3fm 

dyidy2dy3 


dfm  v''  dfm  dm 

=  Z^~nt — >  yi  =  x,t 


l 

Neq 


=E 


diq  Syi’ 

9/m  <9% 


Neq 

E 


92/m  <9?q  <9ufc 


l  dui  dyxdy2  ^  duiduk  dyi  dy2 


(yi,V2)  =  (x,x),  (t,t),  (x,t) 


Neq 

:E 


a/m  53u; 


Neq 

E 


d2f„ 


d2ui  duk  d2ui  duk  d2ui  duk 


l  dui  dy\dy2dy3  ^  duiduk  \dyxdy-2  dy3  dy3dy3  dy2  dy2dy3  dyi 


(7) 


Neq 

+  E 

l,k,p 


d3fm  dui  duk  dup  ^  t  A  _  (x,  x,  x), 


(yi,  2/2, 2/3)  = 


duidukdup  dyi  dy2  dy3  (x,x,t),  (x,t,t), 


for  to  =  1 , ,M.  Equation  Eq.  (7)  shows  the  derivatives  required  by  the  fourth  order  scheme  but  these 
equations  will  continue  to  the  derivatives  required  by  the  desired  order. 

By  substituting  Eqs.  (2)  and  (3)  into  Eq.  (1)  in  a  given  SE (j,n),  we  demand: 


du*m{x,t\j,n)  df^(x,t-,j,n) 


-  Sri 


TO  =  1,  ...  ,  M. 


dt  dx 

Then  by  taking  derivatives  of  Eq.  (8)  in  both  space  and  time  we  get 

«t)j  =  M?  -  (L)",  «*t)"  =  (Sm J”  -  «u)i  =  (Oi  - 

This  result  can  be  written  in  a  more  general  form  as 


(8) 


= 


/* 
rr 


(9) 


for  m  =  1,  . . . ,  M,  z  =  0, . . .  A'm,  *  =  1, . . .  Nm  —  2.  As  shown  in  Eqs.  (7)  and  (9),  the  only  independent 
variables  are  the  spatial  derivatives  for  each  governing  equation.  As  such,  there  are  (Nm  +  1  )M  unknowns 
for  M  equations  associated  with  a  mesh  point.  For  example,  a  fourth-order  representation  of  the  Euler 
equation  would  contain  4  unknowns  per  equation,  giving  a  total  of  12  unknowns  for  the  one-dimensional 
Euler  equations. 

To  proceed,  for  each  of  to  =  1, . . . ,  M,  the  unknowns  are  categorized  as:  (i)  even-order  derivatives, 
,  umx 2n  and  (ii)  odd-order  derivatives,  (umx)j,  n  =  ( OD  -  l)/2. 

In  what  follows,  we  introduce  an  arbitrary  order  CESE  c-r  scheme  for  a  system  of  M  PDEs.  The  c-r  scheme 
uses  space-time  integration  for  the  advancing  formula  for  the  even-order  derivatives,  while  the  odd-order 
derivatives  are  calculated  from  a  central-difference-like  procedure. 
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II. A.  Even-Order  Derivatives 


It  can  be  shown  that  by  differentiating  Eq.  (8)  twice  with  respect  to  x,  we  obtain: 

du*mxx{x,t-,j,n)  ,  df^xx(x,t-,j,n )  d2s*m 

- oJ - H - H - =  -XT’  m  =  l,...,M, 

at  ox  oxz 

or  in  a  more  general  form 

du*mx2z(x,t-,j,n)  df^(x,t;j,n)  d2zs*m  _ 

- di - + - di - "  z  —  0,1, ,  (Nm  —  l)/2. 

Consider  the  Euclidean  space  £2  with  the  coordinates  (^1,^2)  =  (x,t).  Aided  by  defining: 

h *mx2*  d={f^x2z{x,t;j,n),ullx2z(x,t-,j,n))T,  z  =  1, . . . ,  (NM  -  l)/2, 
and  the  divergence  theorem,  Eqs.  (8)  and  (11)  can  be  transformed  into  integral  equations  as: 


Knx2*  '  —  ii2‘  i  2  —  0,  ,  ( Nm  —  l)/2, 


S{V) 

where  S(V)  is  the  closed  boundary  of  an  arbitrary  region  V  and  ds  is  defined  in  Fig.  (2). 


(10) 


(11) 


(12) 


Fig.  2:  Space-time  integration  over  an  arbitrary  closed  domain  V. 


We  define: 


* 

^ mxztk 


r*  def 
J mxztk  ~ 


dz+ku*m 

(  Ax 

dxztk 

V  4 

dZ+kf  *m  , 

(Ax 

dxztk 

{  4  . 

dz+ksm  . 

(  Ax 

dxztk 

K  4  . 

(  Af  \  k 


- 


At' 


U 


(13) 


(14) 


(15) 


where  Ax  =  Xj+ 1/2  —  ^-'j-1/2  and  At  =  tn  —  tn~l .  In  order  to  write  equations  more  compactly,  any  local 
constant  enclosed  within  a  square  bracket  will  be  evaluated  at  the  location  specified  by  the  subscript  and 
superscript  written  on  the  enclosing  square  bracket,  e.g.: 


A  T  A  t. 

i^mxx) j  “1“  i^mxxx) j  ^  H"~  y^mxxt) j 


Ax 


At' 


'U'mxxx  2  ’  ^"mxxt  ^ 


J  3 
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Aided  by  Eqs.  (13)  and  (14),  Eq.  (12)  gives: 


1 


(?w)”  -  Aa. 

Nm  —  z 

)  £ 


Smxz  dV+ 


1  ~  2fc 


2  ^  (fc  +  1)! 


( 

At 

n- 1/2 

i  At 

( 

^jmxkJrZ  T  ^  ^  J  mxztz+k 

+ 

3~  1/2 

(  1)  'UJmxk+z  ^  ~  Jmxztz+k 

n-l/2N 
J+l/2  > 


(16) 


y  _ _ 

it  (2fc  +  i)! 


{umx2k+z)  j 


Equation  (16)  provides  and  explicit  formulation  for  all  even  spatial  derivatives.  As  long  as  the  highest 
even  derivative  is  calculated  first  the  last  term  on  the  RHS  will  have  already  been  calculated.  It  should  also 
be  noted  that  the  *  is  absent  from  the  source  term.  This  is  because  the  source  term  treatment  varies  when 
dealing  with  different  flow  physics. 


II. B.  Odd-Order  Derivatives 

In  order  to  compute  the  odd  derivatives  a  central  differencing  approach  is  applied  following  the  c-t  scheme. 
There  are  two  possible  formulations  for  the  odd-order  derivatives  (i)  the  standard  c-t  scheme  which  is 
applicable  if  there  are  no  discontinuities  present  and  (ii)  a  re-weighted  c-t  scheme  which  is  used  if  there  are 
discontinuities  in  the  flow  field. 

In  order  to  mitigate  the  dissipation  as  the  local  CFL  number  decreases  the  central  differencing  is  applied 
at  points  P+  and  P_.  Where  P±  are  points  located  at 

/\t  /\t 

X  ( P+ )  =Xj  +  (l  +  T)—=  Xj+1/2  ~  (1  -  T)  — ,  (17) 

x  (P~)  =  Xj  -  (1  +  t)^-  =  Xj_1/2  +  (1  -  t)^.  (18) 

Where  r  is  the  absolute  value  of  the  local  CFL  number. 

First  we  define  u^nSz(P±)  to  be  the  Taylor  series  expansion  of  (um2*)ri  from  ( Xj,tn )  to  x(P±).  Then  we 
can  solve  for  umSz+ 1  by  subtracting  u^nSz(P~)  from  u^nSz(P+): 


'UJmxz+1 


rUjmxz  ^mxz{P  ) 

2(1 +r) 


NM~1~z 

2  1 

(2fc  +  i)!Um$2fe+1+z(1  T)  ’ 


(19) 


for  z  =  0, 2, 4, ... ,  Nm  —  1  and  m  =  1,2,. . .  ,M.  Since  we  can  not  calculate  u^Sz  (P±)  we  approximate  it  by 
umi'(^±).  Where  u'mSiz(P±)  is  the  Taylor  series  expansion  from  {Xj±i/2,tn- 1/2)  to  x(P±)  respectively. 


7/ 

^ mx z 


2(1  +t) 


2  1 

5Z  (2k  +  i)iUm$2fc+1+z^  +  r)  fc’ 


When  discontinuities  are  present  in  the  flow  field  a  re-weighting  of  the  derivatives  is  required.  All 
previously  derived  re-weighting  schemes  used  in  second  order  CESE  schemes  should  be  applicable  to  the 
higher  order  CESE  schemes.  In  this  paper  the  derivation  of  the  W2  scheme"  will  be  presented.  First  the 
function  W  is  given  as 


W±{x-,  x+,a) 


k-l“  +  k+lc 


(20) 


To  remain  stable  in  the  presence  of  discontinuities  a  >  1.  The  odd-order  derivatives  are  now  defined  by: 


{V"mxz 


n  def 
3  ~ 


(uJm—  )  2  (S^mx—Z  )  +  (^771+)  2  (^mx+z)  5 


(21) 
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where 


(wm±)z  —  W±(umS_z ,  ,  az),  (22) 

with 

~  def  . 

Mmip  =  ± - —T, - , 

i  +  T 

«mi>T  =f  0"  ^  (Vi-0jTl/2). 

where  (M,mx*-1)j±i/2  Taylor  series  expansion  from  (xj±_i/2,n  —  1/2)  to  (%j±—i/2,n)- 

The  above  equations  provide  an  explicit  formulation  for  the  odd  spatial  derivatives  when  discontinuities 
are  present  in  the  flow  field. 


III.  Jacobian  Matrices 


In  this  section  we  present  the  derived  Jacobian  matrices  for  the  different  flow  physics  used  as  test  cases. 
The  three  flow  physics  used  are  the  Euler,  linear  acoustic,  and  convection  equations. 

To  construct  a  fourth-order  CESE  solver  for  the  one-dimensional  Euler  equations,  we  plug  the  following 
unknown  vector  and  flux  function  vector  into  Eq.  (1): 

U  =  (u!,u2,  u3)t  =  (p,  pv,p/("/  -  1)  +  pv2/ 2)T  , 

F  =  (pv,pv2  +p,  ( pE  +  p)v)T 

=  (u2,  (7  -  1)«3  +  1/2(3  -  7)2x2/221, 7U2M3/M1  -  1/2(7  -  1  )ul/ul)T  , 

where  p  is  the  density,  v  is  the  velocity,  E  is  the  internal  energy  and  7  is  the  ratio  of  specific  heat.  For 
the  higher  order  derivatives  the  order  of  differentiation  does  not  mater,  i.e.  /m,  =  fmk,  1  and  fmi,k,p  = 

=  fmk,ilP  =  fmkiPii  =  fmp,i,k  =  fmp,k ,,  so  only  the  first  such  permutation  is  expressed.  The  first-order 
derivatives  of  the  flux  function  /1  are: 

fh  =  0,  /l2=l,  fU=  0, 


and  the  second-  and  third-order  derivatives  are: 


fh,k=  0,  l,k  =  1,2,3, 
/ii,*,P=0>  l,k,p=  1,2,3. 


The  first-order  derivatives  of  the  flux  function  f2  are: 


hi  -  t. 

(7-3)4  = 

ha  — 

(3-7A 

ha  =  7- 

1, 

2 

“1 

Ml 

the  second-order  derivatives  are: 

hi,i 

_  “2(7-3) 

-  —5?— 

1  hi,'. 

_  (7-3)it2 
>  u\  ’ 

hi,  3  =  0, 

ha, 2 

_  (3-7) 

U\  ’ 

ha,'. 

,  =  0, 

h3,s  =  0= 

and  the  third-order  derivatives  are: 

hi, 1,1  = 

.  0 u2('y~ 3) 

■  — ’ 

> 

_  0(7-3)112 

U1  ; 

>  — 

:0, 

f 2l,2,2  — 

(7-3) 

'  ’ 

f 2l ,2,3 

=  0, 

hi,  3, 3  — 

:0, 

f 22,2,2  — 

:0, 

f 22,2,3 

=  0, 

h3,3,3  — 

:  0. 

The  first-order  derivatives  of  the  flux  function  f3  are: 


f  _  7“2“3  1  (7  — 1)«2  f  _  7“3  3  11  “2  f  __  7“2 

J  3i  -  S — >  J32  -  nr  -  2\1~  J33-—, 
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the  second-order  derivatives  are: 


f  _  niu2u3  o 

J 3i,i  -  - d 

_  o(7-l)“2 

“  - 5 - 5 


,  (7-l)u 


/s2,2  =  -3 

the  third-order  derivatives  are: 


f,  =  -2 
f,  =  -X 

J  j2,3  Uj  : 


o(7-l)t 

7/3 


/3lll=  12^11  -  6^1 


/3l,2,2  =6(7-1)^, 

f  =  07— 1 

J32,2,2  J  51?  > 


/3i,i,2=  2^f-9 

/s  1,2,3  = 
t/*32,2,3 


(7~lM 


/3l,3  =  “ 

/33,3  =  0, 


/3i.i..  =  2^ 

f 3l,3,3  =  0 

h3,3,3  =  0- 


Next  we  present  the  derived  Jacobian  for  the  linear  acoustic  equation 

U  =  (ui,u2)T  =  ( p,v)T  , 


F  =  I  pinfV,  —  P  I  =  I  PinfM2,  -^Ul 


Pinf 


inf  „ 


Pint 


The  first-order  derivatives  for  the  flux  functions  are: 


fli  —  O5/12  —  Pinf 

/2l  =  ^k/2a  =  0. 

Pinf 


Since  all  of  the  first-order  derivatives  are  constant  the  higher  derivatives  are  zero.  This  reduces  the  calculation 
of  all  fluxes  to  a  matrix  vector  multiplication. 

Finally  we  present  the  derivation  needed  for  the  convection  with  source  term  equation. 


du  du  ,  . 

— — b  a—  =  aS 0  cos(a;) 
dt  dx 


where  a  and  Sq  are  constant.  In  this  case  the  source  term  is  only  a  function  of  space  and  an  exact  integration 
is  possible.  For  S  =  aSo  cos(a;)  the  integrals  become 


aS®~Y  sin(a')l*3-i/2 


/2’ 


JJ  Sxs  =  -aSo  ^  Sin(a:)|^+^  , . . . , 

JJ  SS2n  =  (~l)n/2aS0  (J^j  jjsin(x) for  n  =  0,1,... 


Nm  —  1 
2 


IV.  Results 

The  following  test  cases  show  how  the  CESE  method  improves  as  the  order  of  accuracy  of  the  method 
employed  increases.  These  cases  are  presented,  including  (i)  acoustic  waves  modeled  by  the  linearized  Euler 
equations,  (ii)  acoustic  waves  modeled  by  the  nonlinear  Euler  equations,  and  (iii)  a  simple  convection  equation 
with  a  source  term.  For  all  three  cases,  we  calculate  the  order  of  accuracy  by  using  the  following  formula: 


e2 


where  (f>i  is  defined  as  the  difference  between  the  analytical  and  numerical  solution  and  Ax,  is  the  grid 
spacing  at  a  given  location  i.  In  all  cases  Ax  is  constant.  The  rate  of  convergence  is  taken  as  the  slope  of 
the  best  fit  line  through  the  points  (logio(Ax),  logi0(£2)  ). 


7  of  14 


DISTRIBUTION  A. 
Approved  for  public  release; 
distribution  unlimited. 


American  Institute  of  Aeronautics  and  Astronautics 


The  first  test  case  is  the  calculation  of  acoustic  waves  by  solving  the  linearized  Euler  equations.  The 
analytical  solution  for  the  linearized  acoustic  wave  equation  is 

epinf  /2ri7r  \ 

p  =  Pinf  H - cos  —^{x  -  ainf t) 

&inf  \  l  J 

(  2  mr  \ 

U  =  c/inf  +  £COS  I  — —  (X  -  dinft)  I 

„  l  l 

for  -  <  x  <  - ;  t  >  0 

2  2 

where  p,  U ,  a,  n,  l,  and  e  are  respectively  the  density,  velocity,  speed  of  sound,  number  of  waves,  length  of 
the  domain,  and  an  amplification  factor.  The  speed  of  sound  is  equal  to  y/j p/p  with  7=1.4.  The  values  with 
a  subscript  inf  are  mean  values  of  the  flow  variables.  For  this  test  case  p;nf=l,  Pinf=7)  £=10  ,  n=  1,  and 

1=2.  The  run  time  was  equal  to  4.25^2--  which  allows  the  wave  to  propagate  through  the  domain  4.25  times. 
The  CFL  number  is  constant  throughout  the  domain  and  is  equal  to  0.75.  As  seen  in  Fig.  (3)  our  desired 
order  of  convergence  closely  matches  the  actual  order  of  convergence.  Table  (1)  shows  the  desired  order 
of  convergence,  the  actual  order  of  convergence,  and  the  normalized  time.  The  relative  numerical  cost  was 
calculated  by  taking  the  average  simulation  time  per  cell  per  iteration  for  multiple  resolutions  and  dividing 
it  by  the  cost  of  the  2nd  order  version.  The  computational  scaling  shows  that  by  doubling  the  order  of  the 
Taylor  series  the  time  required  to  complete  a  simulation  will  increase  by  about  22  2. 


logio  (  A  x) 

2nd  Order  -  12*^  Order  - 

4th  Order  -  16th  Order  - 

8th  Order  -  20^  Order 

Fig.  3:  The  £2  norm  of  numerical  solutions  where  points  are  actual  calculated  data  and  the  line  is  a  best-fit 
cure  of  the  data. 


Desired  Order 

Actual  Order 

Normalized  Time 

2 

2.03 

1.00 

4 

4.06 

3.43 

8 

8.12 

14.47 

12 

12.21 

34.33 

16 

16.48 

67.98 

20 

20.52 

116.89 

Table  1:  Convergence  rates  for  the  linear  acoustic  solver  and  the  average  normalized  time  for  case. 


The  above  cases  showed  that  we  were  able  to  achieve  higher-order  convergence  for  coupled,  linear,  wave 
equations.  In  the  following,  we  show  fourth-order  convergence  for  solving  non-linear  equations  but  for  linear 
physics  .  The  same  analytical  solution  used  for  the  linear  acoustic  equation  is  used.  There  are  some  problems 
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when  using  this  solution  because  the  Euler  solver  will  capture  the  non-linearities  present  in  the  flow  field 
while  the  “analytical”  solution  does  not.  This  will  lead  to  increasing  errors  in  the  analytical  solution  as 
A  a:  decreases.  To  mitigate  the  error,  the  perturbation  was  reduced  to  10~6.  For  this  test  case  Pi„f=l/7, 
Pinf=l,  n= 2,  and  1= 4  and  the  simulation  time  is  2.5-^.  The  CFL  number  is  almost  constant  throughout 
the  domain  and  is  equivalent  to  0.8.  As  seen  in  Fig.  (4)  we  achieved  a  convergence  rate  of  2.01  and  3.97  for 
the  second-  and  fourth-order  CESE  schemes,  respectively. 


Ath  Order  norm  +  2nd  Order  l2  norm  + 

Ath  Order  -  2nd  - 


Fig.  4:  The  l2  norm  of  the  numerical  solutions  of  the  Euler  solver  for  solving  the  acoustic  waves. 


The  final  test  case  is  solving  the  convection  equation  with  a  source  term.  Under  the  periodic  boundary 
condition,  the  analytical  solution  to  this  problem  is 

u(x,  t)  =  cos{x  —  at)  +  S'osm(x), 

— 27T  <  x  <  27 r;  t  >  0. 

For  the  convergence  tests,  a  =  Sq  =  1  and  the  test  time  was  set  to  2.5-,  where  l  is  the  length  of  the  domain. 
In  all  calculations,  CFL  =  0.7.  Shown  in  Fig.  (5)  and  Table  (2),  the  actual  convergence  rate  agrees  well  with 
the  order  of  accuracy  of  the  scheme  employed.  The  computational  scaling  shows  that  by  doubling  the  order 
of  the  Taylor  series  the  time  required  to  complete  a  simulation  will  increase  by  about  22'2. 


Desired  Order 

Actual  Order 

Normalized  Time 

2 

2.00 

1.00 

4 

4.01 

5.22 

8 

7.95 

23.65 

12 

12.12 

59.04 

16 

16.05 

115.83 

20 

20.09 

190.95 

Table  2:  Convergence  rates  of  the  numerical  solutions  of  the  convection  equation,  and  the  averaged,  normal¬ 
ized  time  for  case  with  different  order  of  accuracy. 
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2th  Order  -  12th  Order 

4th  Order  -  16th  Order 

8th  Order  -  20th  Order 


Fig.  5:  The  £2  norm  of  numerical  solutions  of  the  Convection  equation  with  source  term.  The  symbols 
represent  the  actual  calculated  data  and  the  lines  represent  the  best-fit  curves  of  the  data. 


Another  important  aspect  to  consider  is  whether  the  higher-order  resolution  is  worthy  for  the  additional 
computational  cost.  For  this,  we  refer  to  Figures  (6)  and  (7).  Shown  in  the  figures,  it  is  more  efficient  to  use 
a  higher-order  method  rather  than  increasing  the  resolution  for  a  linear  solver. 


CPU  time  (seconds) 


2nd  Order  * — 12th  Order  - 1 — 

4th  Order  - ' -  16th  Order  - ' - 

8th  Order  - ' -  20th  Order 

Fig.  6:  The  £2  norm  versus  the  computational  time  for  the  solutions  of  the  linear  acoustic  wave  equations. 
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2nd  Order  - 1 - 

12th  Order 

4t,!  Order 

16th  Order 

8th  Order  - ' - 

20th  Order 

Fig.  7:  The  f2  norm  versus  the  computational  time  for  the  solution  of  the  convection  equation. 


Next,  we  demonstrate  the  new  high-order  CESE  method  by  examining  numerical  resolution  for  shocks 
and  contact  discontinuities.  We  will  run  two  different  test  cases  at  various  resolutions  and  compare  the 
results  with  an  “analytical”  results  obtained  by  using  very  fine  mesh.  Three  solvers  will  be  used  in  the 
test  cases,  the  second-  and  the  fourth-order  CESE  and  the  fifth-  order  space  third-order  time  monotonicity 
preserving  (MP53)  method. 1  The  comparison  between  the  4th-order  CESE  scheme  and  the  MP53  method 
is  not  completely  valid  because  the  CESE  scheme  employed  has  higher  order  in  time  but  low  order  in  space. 
The  test  cases  chosen  are  Woodward’s  blast  wave  problem  and  Shu  and  Osher’s  problem.  Woodward’s 
blast  wave  problem  consists  of  two  shock  waves  of  different  strengths  heading  towards  each  other  with  wall 
boundary  conditions.  The  initial  conditions  are 


(u,P,p)  = 


(0, 

1, 

103) 

0  <  x  <  0.1 

(0, 

1, 

10-2) 

0.1  <  x  <  0.9 

(0, 

1, 

102) 

0.9  <  x  <  1.0 

0.0  <  t  <  0.038. 

The  second  test  case  is  Shu  and  Osher’s  problem,6  in  which  a  Mach  3  shock  moves  to  the  right  and  collides 
with  an  entropy  disturbance  moving  to  the  left.  The  boundary  conditions  are  non-reflective  and  the  initial 
conditions  are 


(u,P,p) 


(2.629369,  3.857143,  10.3333) 

(0,  1  +  0.2  sin(57rx),  1) 


—  1  <  x  <  —0.8 
—0.8  <  x  <  1.0 


0.0  <  t  <  0.47. 
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0.5  0.55  0.6  0.65 


0.7  0.75  0.8  0.85  0.9 

x 


MP53  >  CESE  2nd 


CESE  4th  Converged 


Fig.  8:  Plots  of  the  density  profiles  of  the  Woodwards  blast  wave  problem.  The  converged  simulation  was 
done  by  using  the  fourth-order  CESE  with  3201  points.  For  better  presentation,  only  a  subset  of  the  domain 
is  shown  in  these  plots. 
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X 


MP53  «  CESE  2nd  *  CESE  4th  Converged 

Fig.  9:  Plots  of  the  calculated  density  profiles  for  Shu  and  Osher’s  problem.  Each  plot  has  a  different  spatial 
resolution.  The  converged  simulation  was  done  by  using  the  fourth-order  CESE  scheme  with  3201  points. 


For  both  simulations  using  the  second-  and  fourth-order  CESE  methods,  a  in  Eq.  (20)  is  set  to  be  unity. 
Shown  in  Figures  (8)  and  (9),  the  fourth-order  CESE  method  provides  more  accurate  solutions  than  the 
second-order  CESE  method  in  the  region  where  the  solution  is  smooth.  When  discontinuities  are  present, 
the  fourth-order  scheme  still  does  a  better  job  but  has  more  overshoots  than  the  second-order  CESE  method. 
Moreover,  the  results  obtained  by  using  CESE  schemes  compare  favorably  with  that  by  the  MP53  scheme. 

V.  Concluding  Remarks 

In  this  paper,  we  extended  Chang’s  fourth-order  CESE  method  for  one  convection  equation  for  solving 
a  system  of  coupled  hyperbolic  PDE’s  with  arbitrarily  high-order  convergence.  Numerical  results  show  that 
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the  extended  algorithm  can  achieve  higher-order  convergence  for  both  linear  and  non-linear  hyperbolic  PDEs. 
The  shock-capturing  capability  of  the  new  method  was  comparable  to  that  of  the  original  second-order  CESE 
scheme  as  well  as  that  of  the  fifth-order  space  third-order  time  monotonicity  preserving  scheme.  Further 
development  of  the  high-order  CESE  method  can  benefited  from  investigations  in  several  areas,  including 
the  effect  of  different  limiters  on  the  higher-order  derivatives,  and  the  effects  of  the  boundary  condition 
treatments  and  the  source-term  treatments  for  high-order  accuracy. 
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